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Abstract 

Using the prescription of the nonequilibrium statistical operator method, we derive a non- 
Markovian generalization to the kinetic theory described by Walser et ai. [Phys. Rev. A 59, 
3878 (1999)]. Quasi-particle damping and effects arising from the finite duration of a collision 
are introduced to include terms beyond the Born approximation. Such a self-consistent theory is 
shown to conserve energy to second order in the interaction strength, even in the Markov limit. 
This kinetic theory is applied to a simple model of a Bose gas confined in a spherical trap to study 
the full real-time evolution towards equilibrium. A modified form for the damping function, is seen 
to strongly improve the energy conservation. Based on a linear response calculation, we predict 
the damping rates and frequencies of the collective excitations. We demonstrate the emergence of 
differing time scales for damping and equilibration. 
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I. INTRODUCTION 



In the mean-field approximation, a Bose-condensed phase is well described by the Gross- 
Pitaevskii (GP) equation Examples of collective effects that arise at the mean- field 
level include the formation of vortex states [0, ^ ||, |^ and collective excitations |]^, ||, 
0. The GP equation is often sufficient to describe the dynamics at T = 0. However, 
at finite temperatures, T near [the Bose-Einstein condensation (BEG) temperature], 
effects due to the presence of a thermal component may be important. The hot thermal 
atoms interact with each other and with the condensed atoms via binary collisions. In 
fact, coUisional dynamics is the essential mechanism for evaporative cooling leading to BEG. 
Finite temperature effects generate phenomena such as phase diffusion and damping of 
collective excitations. In order to describe these effects, inclusion of collisions arising from 
the thermal component becomes necessary. A generalized GP equation, for example the 
Hartree-Fock-Bogoliubov approach, includes collisions only as first-order energy shifts and 
is therefore valid only at low temperature. In the other limit, T > Tc (high temperature), 
the condensate component may be negligible and the dynamics is then dominated by the 
thermal component. This limit is well described by the quantum-Boltzmann (QB) equation. 
A second order kinetic theory, which reduces to the GP and the QB equations in the limits 
T = and T > T^, respectively, is necessary to provide a complete description. 

Several different versions of quantum kinetic theories have been formulated to describe a 
quantum degenerate Bose system [^, 0, |T3|, |T5|, . However certain questions 

remain unanswered. An important issue is the question of Markovian versus non-Markovian 



dynamics in a non-homogeneous system. Previous attempts to address this issue |T^, [1^, 
pO| , plj p^] , have been limited to the discussion of system with translational symmetry. 

In the case of a sufficiently dilute nonhomogeneous system, one may consider the neglec- 
tion of memory effects under the principal of the rapid attenuation of correlations. Such an 
ansatz will result in a Born-Markov approximation. While in this limit kinetic equations 
are vastly simplified, they do not conserve energy. Also, the Born-Markov theory involves a 
regularizing factor 77 — [^S], |2^. Such a limiting procedure is computationally nontrivial 
to implement numerically if one is dealing with a finite system with discrete energy levels. 
If the resolution is instead to include memory effects and initial correlations, then the full 
non-Markovian theory which arises is described by complex kinetic equations. 
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In this article we wish to address some of these issues by generahzing the kinetic theory 



in Ref. |]T0| to include effects beyond the Born approximation such as quasi-particle damping 
and effects arising from the finite duration of a collision. We hypothesize the form of the 
damping function by looking at the corresponding retarded Green's function pB| , p7[ . We 



show that such a perturbation theory conserves energy up to second order in the interaction 
strength. 

We apply this theory to a simple model of a Bose gas confined in a spherical box above 
Tc- The simplicity of this model allows us to study the real-time (nonequilibrium) evolution 
of the system from some initial state to a final state of equilibrium. The equilibrium state 
is essentially a self-consistent steady-state solution of the kinetic equation. We perform a 
linear response calculation in order to verify the stability and to study the damping rates 
and frequencies of the steady-state solution. We show the existence of different time scales 
corresponding to damping and equilibration. 

This article is organized as follows. In Sec. II we derive a non-Markovian generalization 
to the kinetic theory of Ref. [|10] using the prescription of a nonequilibrium statistical op- 
erator method . Both quasi-particle damping and damping arising from the finite time 



of collision events will be discussed in Sec. |ITT[ This has implications on the underlying 
symmetries of the theory and their associated conserved quantities. In Sec. |V|, we apply 
this generalized kinetic theory to a simple model of a dilute Bose gas confined in a spherical 
box and obtain a self consistent steady-state solution to the second-order kinetic theory. 
Finally in Sec. 0, we study the response of the system to a small perturbation. This will 
help us to predict frequencies and damping rates of collective modes. 

II. KINETIC EQUATIONS 

We start with the many-body Hamiltonian for a weakly interacting Bose gas given by 

H = H^'^ + H^'\ (1) 
where if'-^-* is the single particle Hamiltonian that is defined as 

H^'^=e'^ala, (2) 



3 



using the implicit summation convention for repeated indices. The two-body energy IS 
given by 

^(1) ^ 01234 ^{0^0304. (3) 

The bosonic operators di and d\ annihilate and create a particle in a single particle state 
|1), respectively. The abbreviated notation |1) represents a state specified by a complete set 
of quantum numbers for both the motional and electronic degrees of freedom. 

We assume that the particles are confined by an external trapping potential l^xt- Thus, 
the matrix elements of the single particle Hamiltonian are given by 

e'' = (l||^ + Kxt(i)|2), (4) 

where m is the mass. The binary interaction is mediated by a short range repulsive potential 
Vbin- This gives the symmetrized (S) matrix elements 



il234 



® (2|Vbin(x® 1 - 1 ® x)|3) ® |4). (5) 



In the low energy limit Vbin can be approximated by a contact potential with the matrix 
elements given by 



^1234 _ 

2 



/oo 
(l|x)(2|x)(x|3)(x|4)d3a;. (6) 
-00 



The interaction strength Vq, is related to the scattering length as by Vq = A-nfi^aa^/m. 

Here, we use the well-known nonequilibrium statistical operator method ^ to ob- 
tain an explicitly non-Markovian version of the kinetic theory |T^. In this approach, the 



nonequilibrium state of a weakly interacting quantum gas is specified by a set of single time 
master variables. For our system the most important master variable is the single particle 
density matrix /(t), 

/i2(t) = (4ai) = Tr \ a\a^ a{t) \ , (7) 



where (...) = Tr{... (T(t)} and a{t) is the statistical many-body density operator. We 
focus our studies on the temperature regime above and in close proximity to the critical 
temperature for BEC. For this reason, we do not consider either symmetry breaking fields, 
(ai), or the anomalous fiuctuations, (0-102). We therefore define = l,7jr. = d\ai\k E 
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{(1, 2)}} as our complete set of relevant operators. The expectation values of these operators, 
lk{t) = ilk), are the only quantities that will appear in the final kinetic equations. 

The time evolution of the nonequilibrium statistical operator a{t) is described by the 
Liouville equation with an extra source term on the right hand side: 

Jt^it) + ^[H,Ht)] = -V (8) 

Such a source term breaks the time reversal symmetry of the Liouville equation and rep- 
resents a convenient way to incorporate the irreversible character of macroscopic processes. 
The relevant distribution a^^\t) given by 

^^'\t)=a^;!;^,^^=exp{%T\t)}, (9) 

where T'^(t) are the Lagrange multipliers, represents a special solution that maximizes the 
information entropy S' = — Trjcr' log(cr')} for the given averages Jkit)- Furthermore, at 
some initial instance t = to in the remote past, we can assume that a(to) corresponds to its 
noninteracting value and therefore 

a{to) = a^'\to). (10) 

The Lagrange multipliers T^{t) are calculated from the self-consistency condition 

7fc(t) = Tr{7fc(T(°)(t)} = TT{%a{t)}. (11) 

This essentially enforces the Chapman- Enskog condition [^] for the restricted set of relevant 
operators at all times. 

From the Liouville equation (P) one can easily establish the basic equations of motion for 
the average values, 

+ ^Tl{|7/<",%]a(()}. (12) 

The form of enables us to express the first trace on the right hand side of the above 
equation in terms of the averages Jkit)- The second trace plays the role of the "collision" 
term, evaluation of which requires us to seek an integral solution of the Liouville equation. 
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But, before we proceed, it is instructive to repartition the total Hamiltonian Eq. (|l|) into 
single particle and two-particle contributions. 



^ (0) ^ (1) 

H = H it) + H (t) 



H^'^ + Qit)] + \h^'^ - Qit) 



(13) 



This modification anticipates self-energy shifts 



Q{t) = Q''{t)ala,, 



(14) 



which will inevitably arise in the course of the calculation. 

An integral solution for a{t) can then be obtained easily from the Liouville equation 
(1), by using the single particle time evolution operator ?7^°''(t, ti), the boundary condition 
Eq. (^, and an additional partial integration. Thus, one finds 

a{t) = a(°)(t)- f dhe-'^^'-''^U^'\t,t,) 

J to 



X 



(15) 



.(0), 



Since the Hamiltonian H (t) depends on time through the Q(t), the time evolution operator 

1/(0) 

is in general a time-ordered exponent 



f/(o)(t,to) = Texp 



h 



^ (0) 

dtiH [ti] 



to 



(16) 



To establish the time derivative da^^\ti)/dti in Eq. (|T5D , we recall from Eq. that the 
relevant operator depends only implicitly on time through the averages 7/c(ti). Moreover 
by exploiting the transformation properties of a quantum Gaussian (see Appendix one 
finds that 



+ lTr|[F\t),7.]a(t)U,,aW(t). 



(17) 



Using these relations we eliminate the time derivative of a^'^^ti) in Eq. (|T5|), to obtain the 
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integral form of the statistical operator 



^ J to 



X 



- (1) , 



(18) 



Since we are only interested in a weakly interacting gas, we seek a power series expansion 
in the interaction strength 



(19) 



a 



(t) = -1 / dtie-''(*-*^)f/(°)(t,ti) 



to 



X 



- (1) , 



(20) 



With this explicit expression for the statistical operator, the evaluation of the equation of 
motion (O) is straight forward and one obtains the quantum-Boltzmann equation 



|/(t) = C[f] + C[fl 



(21) 



Here, Wick's theorem has been used to express the higher order averages in terms of the 
single particle ones. The kinetic operator C consists of a reversible Hartree-Fock (HF) part 
£hf and a collisional quantum-Boltzmann contribution Cqb, 



m 



CnF[f]+CQM, 



HHF{t)fit), 



- r(i+/)(i+/)//- 



-f^HF = £ + 2 Uf is the Hartree-Fock Hamiltonian, with f/j^ = 2 0^^^^/32 the self energy, while 
F's are the collision integrals given by 



^ABCD ~ ^ J '^^1 



-»?{t-ti)^1234 ,1"2"3"4" 



X /Ci"i'(t, tl)A^2"2'(^5 ti)A^3/;3/(t, ti)IC\„^,(t, ti) 
X v431'(tl)542'(il)C4'2(tl)/^3'5(tl) 



(22) 
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The propagators Kit, to) are given by 

Kit, to) = Texp 



dti HuF{ti 



to 



(23) 



Therefore, unhke the colhsion terms of Ref. ||T0[, the F's defined above depend on the past 
history of the system. Thus Eq. (p^) represents a non-Markovian generahzation of the kinetic 



equation previously derived in Ref. |jTO|. 



III. CONSERVATION LAWS AND QUASI-PARTICLE DAMPING 

The conserved quantities for a closed isolated system are the total energy E and the 
total number A^. These quantities then represent the constants of motion for the full kinetic 
equation. The total number operator can be represented as a linear combination of the 
relevant operators, 

N = d\a^, (24) 
and therefore the functional N{f), representing the total number is given by 

N{f) = Tr{/}. (25) 
The kinetic equation for A^(/) can then be written as 

|iV(/) = TT{C[f]+C[f]^} 

= Tl {CHF[f]+CHF[f]^} 

+ TT{CQB[f]+CQB[f]^}. (26) 

The first order term on the right-hand side involves the trace of a commutator and is trivially 
equal to zero. The F's associated with the second order terms have the following property 
due to the symmetries of 0, 

Tr {r//{i+/)(i+/) } = Tr {r(i+/)(i+/)//}* . (27) 

As a result, the second order contribution in Eq. (|26|) can also be shown to be zero. Hence 
the total number is a constant of motion. The important point here is that the total 
number conservation is a result of the symmetries of (p and does not dependent on the non- 
Markovian nature of the collision integral. Therefore a Markov approximation would leave 
this conservation law unchanged. 
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While the number conservation is a natural consequence of the self-consistency condition, 
the total energy conservation is not obvious as the Hamiltonian H cannot be represented 
as a linear combination of the relevant operators. We start with writing the total energy 
functional E{f) as a perturbative expansion in 0, 

E{f) = Tr{if(°V(°)} 

+ [Tr{^WaW} + Tr{i^Wa«}" 

+ [Tr{i7(iV(^)} + Tr{i^(°V(2)}j + ... (28) 

Using the self-consistency condition Eq. (0) for the relevant operators, one can show that 
the third and the fifth trace terms in the right-hand side of the above expression drop out. 
The kinetic equation for the energy functional E{f) can then be written as 
d 



dt 



E{f) = TT\ef + Ujf + Ujf 



(29) 



d 



r 



Again, we use the symmetry properties of (p, to write the simplified equation 



d_ 
di 



E{f) 



■ V 



{r//(i+/)(i+/) - r(i+/)(i+/)//} . 



(30) 



One can now see that the energy is conserved only in the ?7 — limit. In this limit the 
kinetic equation (^) represents the Born approximation. A finite value of t] could then be 
thought of as resulting in additional terms that are beyond the Born approximation. Such 
terms model the duration of collision effects and quasi-particle damping. In principle if this 
effect is treated self-consistently, rj will be a time dependent function at least of order (j). 
This means that the rate of change of E given by Eq. ( pO| ) is of the order 0^. 

Also, if one is only interested in times grater than the correlation time Tcor, the finiteness 
of 7] allows us to extend the lower limit of the collision integral to — oo. Now we can 
approximate the f(t') in the non-Markovian expression of T by its instantaneous value f{t) 
to obtain the Markov form 



■.{m)i5 
ABCD 



dti 



-'7(*-*i)^1234 ,1"2"3"4" 



X /Ci"i'(t, ti)}C2"2'(t, ti)/C3„3,(t, ti)ICl„^/(t, ti) 

X Asr{t)Bi2'it)Ci'2it)Dy5it) 



(31) 
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One can verify that the above Markov form results in an energy conservation, up to the 
most significant order given by 



d_ 

di 



E{f) = Tr {/7hf + } . (32) 



If we compare this with the expression for the correlation energy Ecor in Refs. |T^, [19 
the right hand side of Eq. (0) is exactly —dEcor/dt. This is not surprising and can be 
understood more intuitively by writing the collision integral as a sum of two contributions: 
correlation and collision 

r = Teor + r,oi = [ ...dt' + [ ...dt'. (33) 

^-oo Jo 

For a finite r], the Fcor contribution decays to zero as e~^*. Therefore the decaying correlation 
energy Ecor associated with this part shows up in the rate of change of the total energy E. 
In the Born approximation Fcor is constant because 1] —>■ 0. 

The exponential damping results in a widened delta function and therefore the rate of 
change of E(f) Eq. ( p2D can be shown to be of order r]T^"^\ Thus, by including terms beyond 
the Born approximation, we have obtained a collision integral that is Markovian and still 
conserves energy up to 0^ order. Now if we assume that the equilibration time is of the 
order l/F^™^ , the total change in energy AE is therefore 



AE = E{r^) - E{f 
' 1 

p(m) 

where /™ and f^"" are the initial and equilibrium distributions, respectively. 

The importance of the damping term becomes more obvious when one attempts to solve 
the kinetic equation numerically. One no longer has to worry about the 77 limit in the 
Born-Markov approximation, which for a finite system with discrete levels will result only 
in exchange collisions and hence no equilibration. 



IV. APPLICATION TO A DILUTE BOSE GAS IN A SPHERICAL BOX TRAP 

In the previous section, we introduced the general methods and concepts to describe a 
weakly interacting Bose gas under nonequilibrium conditions. We will now apply these to 
a simple model of a typical ^''Rb experiment, as realized by many laboratories around the 
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FIG. 1: Hartree-Fock energies as a function of particles number in the box ground state, /n with 
all other /jj's equal to zero. 



world, for example Refs. J^, The physical parameters are usually quoted in the 

natural units for a harmonic oscillator trapping potential, i.e., the angular frequency uj = 
2tt 200 Hz, the atomic mass mgr = 86.9092 amu, the ground state size cih = [^/ i^sr)]^^'^ = 
763 nm, and the s-wave scattering length as = 5.82 nm. 

However, in the present article, we do not pursue the usual harmonic confinement, but 
rather explore the properties of a radial box as a particle trap. This choice is motivated 
by previous studies of the self-consistent Hartree-Fock single-particle states As soon 
as repulsive mean-field potentials are added to the bare harmonic trapping potentials, the 
corresponding eigenstates widen in size and look remarkably close to the eigenstates of a 
box, provided the spatial extensions of the box is chosen appropriately. In particular, we 
pick a box of radius R = 1000 as = 5.82 fim. 

Thus, our model is represented by a spherical trap with box potential given by 



v;xt(r) 




(35) 



The eigenfunctions are a product of spherical Bessel functions and spherical harmonics 



(Im)- 



1pinlm}{r,9,(t)) 



r> R. 



(36) 
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A normalization constant is represented by M[ni)- The eigenenergies are given in terms of 
the wave- vectors which can be obtained from the n-th nodes of the spherical Bessel- 
functions of angular momentum /, 

= {Rk^^i)/Txf Eq. (37) 

Here, Sq represents the ground state energy and defines the energy scale of the problem 

All the physical parameters are scaled with respect to this energy unit and the radius of 
the box R. 

For simplicity, we assume all the atoms to be in the / = state initially. One therefore 
needs to consider only the / = manifold and get for the normalization constant A/(no) = 
n(7r/(2i?^))^/^ and radial wave-vector i?A;(„o) = nn. If the cloud is relatively cold, then most 
of the population resides in the lowest few energy states. Therefore, we can also limit the 
number of radial modes 1 < n < rzmax- For the present case, we take nmax = 5. Obviously, 
all these simplification reduce the number of degrees of freedom significantly and thus we 
are able to study certain aspects of the nonequihbrium dynamics of the trapped Bose gas 
in great detail. 

With the above definitions, the bare single particle box Hamiltonian is given by 

= 8pq (no implicit sum over q). (39) 

The interaction part H'^^^ involves the matrix elements of the interaction potential defined 
by Eq. (I) 



pqrs 



— ^ I sin (px) sin (gx) sin (r x) sin (sx)^, 
vr Jo 

(40) 

which in general have to be computed numerically. Interestingly, in the case of a spherical 
box this integral can be evaluated analytically and simplified to a finite sum of sine integrals 
and cosine functions (see Appendix ^ for details). 

The nonequihbrium state of the above system is represented by the single particle dis- 
tribution function / with time dependence given by Eq. (^). Neglecting the second order 
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collision terms, the first order evolution is governed by the Hartree-Fock Hamiltonian H^p 
and is given by 

p = CHF[f]+CHF[f]^. (41) 

The energy eigenstates of the interacting system are therefore shifted from the bare box 
states due to the self energy effect. These shifts can be significant depending on the total 
particle number. This is clear from Fig. ^ where we plot the eigenenergies as a function of 
/ii, the total particle number in the box ground state (with all other fi/s equal to zero). 

Now note that the time dependent Hartree-Fock equation (|4l|) for the density matrix, /, 
is nonlinear and hence we seek a self consistent solution such that 

f = J2Pi^^M{e^\, (42) 

i 

where i^Hpl^i) = and for a Bose-Einstein distribution P{e) is given by 

= — rf r (43) 

exp((£: - nj/kBT) - 1 

For a given total particle number, A^, and temperature, (3 = l/ksT, a self-consistent chem- 
ical potential, fi{P,N), and hence a self-consistent Bose-Einstein distribution is obtained. 
For example, let us consider three different total particle numbers = 10, 100, 500 at two 
different temperatures, /? = l/ksT = 0.01,0.5, corresponding to hot and cold clouds, re- 
spectively. In Figs. H and |], we plot the self-consistent solution and the self-energy density in 
the position space representation. We see that the self-energy density is proportional to the 
number density only near the center of the trap and drops off faster with increasing radius. 
This can be attributed to the restricted number of basis states used in our calculation and 
effectively gives a finite range to the two body potential Vbin- 

Up to first order, the / equation is totally reversible. The inclusion of the second order 
terms (collisions) break the reversibility, and therefore represents a relaxation of the system 
from some initial state to a final equilibrium state. Here we will be using the Markov form 



(pl|) for the collision integral. We take the Hartree-Fock self-consistent state for the initial 
condition. 

In the previous section, we interpreted the function exp(— r^r) to account for duration 
of collision effects and quasi-particle damping. But the exponential form was originally in- 
troduced to break the time reversal symmetry and it has the correct long time behavior. 
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FIG. 2: The mean-field density, Uf{r,r), as a function of the radial distance at a relatively hot 
temperature, (3 = 0.01 is shown for three different values of particle numbers: iV = 10 (solid curve), 
N = 100 (dashed), and N = 500 (dot-dashed). Inset shows the corresponding number density as 
a function of radial distance. 
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FIG. 3: The mean-field density, Uf{r,r), as a function of the radial distance at a relatively cold 
temperature, P = 0.5 is shown for three different values of particle numbers N = 10 (solid curve), 
N = 100 (dashed), and N = 500 (dot-dashed). Inset shows the corresponding number density as 
a function of radial distance. 



However, an exponential damping will result in a Lorenzian line shape for the final equilib- 
rium distribution. Due to the long-reaching wings of the Lorenzian curve, in the Markov 
limit, off-the-energy shell collisions get weighted strongly. To seek an improved damping 
function that will have the correct short and long time behavior, we use the equivalence of 
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kinetic theories based on the Green's function approach |^ and the non-equihbrium sta- 
tistical operator method as shown in Ref. The behavior of the retarded Green's 
function g{t,ti) for very large and very small time scales is given by 



g{t,ti 



e-vit-h) t-ti:^T,,, 



(44) 



Therefore, behavior over the intermediate time scale will be best represented by an inter- 
polating function. This is also true for the damping function. From Fig. ^ we see that the 
function 



^(t,ti) = l/cosh(r/(t-ti)), 



(45) 



has exactly this behavior and therefore represents a better choice than the exponential form. 

With either choice of the damping function, for a particular value of the parameter 77, 
a time propagation results in a self-consistent steady-state solution, f'^'^. Figure ^ shows 
such a time evolution for = 500 particles with an initial temperature corresponding to 
(3 = 0.01. We have chose t] to be of the order Re[r] ^ 2.3. As mentioned earlier we see 
that the exponential damping function results in significant transfer of population to the 
excited states. This effect is less with the hyperbolic secant damping function. We get 
different steady-state solutions because different damping functions correspond to different 
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FIG. 4: Comparison between different damping functions, exp(— j^r) (solid), exp(— ryr^) (dot- 
dashed), and 1/ cosh(7yr) (dashed). Note that the hyperbolic secant function asymptotes to an 
exponential form for large r and a Gaussian form for small r. 
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FIG. 5: Evolution of / toward a self-consistent steady-state solution starting from an approximate 
solution. Two different damping functions are used in the evaluation of the F's. Dashed and solid 
curves correspond to exponential and hyperbolic secant respectively. 
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FIG. 6: Linear behavior of 9{fnn) as a function of s^^ ■ Initial distribution shown with diamonds. 
Final distribution shown with circles for the case of exponential damping function, and squares for 
the case of hyperbolic secant damping function. 

initial correlations. Also a plot of 9{fnn) = ln(l//ran + 1) vs the Hartree-Fock energies ef^^, 
shown in Fig. ^ shows that the f^^ is very close to a Bose-Einstein distribution. The slope, 
which represents the self-consistent value of (3, shows that the change in temperature is 
far greater for the exponential damping function compared to the hyperbolic secant case, 
because off-the-energy shell effects are larger as explained previously. 

Thus we have obtained a self-consistent steady-state solution to the second order kinetic 



0.2 
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equation. We emphasize here that the steady-state solution is a result of the real-time non- 
equilibrium evolution of the system. Deriving such an equilibrium solution, whose absolute 
value is plotted in Fig. 0, is a prerequisite step to the study of collective modes and damping 
rates of a dilute gas. 

V. REAL TIME RESPONSE TO PERTURBATION 

The properties of the equilibrium solution f^"^ exhibit the expected characteristics of a 
Bose-Einstein distribution. In order to verify the stability of this solution and to study the 
damping rates of the collective excitations, we will now examine the real-time response of the 
system to a perturbation. First, we will outline the linear response theory and discuss the 
structure of the modes, their frequencies, and the life-time of the excitations. Subsequently, 
we will use these modes to initially prepare the system and to evolve the full nonlinear 
quantum kinetic equation towards equilibrium. 

One of the fundamental properties of the quantum kinetic equation (pT]) 

|/(t) = C[f] + C[f]\ (46) 

is its Hermitian structure. Thus, if we prepare a physical state initially, it will remain 
Hermitian with f{t) = /(t)^ indefinitely. We will now consider a weak perturbation of an 
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FIG. 8: Non-negative frequency eigenvalues shown with crosses for the Hartree-Fock equation and 
squares for the quantum-Boltzmann equation. The dotted hues correspond to difference energies 



sf^- The modes labelled (a) and (b) - non-zero frequency and damping rate, (c) - zero 



mode, and (d) - zero frequency and non-zero damping rate will be considered for further discussion. 



equilibrium state, 



fit) = r + sf{t), 



(47) 



and calculate the first order response of the system. As usual, we want to assume that we can 
decompose a general perturbation into fundamental damped and/or oscillatory eigenmodes 
of the system. Therefore, such a specific perturbation can be parameterized as 



6 fit) = e-^-*5/W + h.c., 



(48) 
(49) 



where denotes a positive frequency amplitude. It turns out to be useful to decompose 
it into quadrature components 



sf!^ 
5f!^ 



1 

2 
1 

2i 



sf^-"^' + sfL-"^ 
sf^-"^^ - Sfl^^ 



(50) 
(51) 



and evaluate the kinetic operator C[f] only for such Hermitian arguments. Prom a Taylor 
series expansion of the kinetic operator around the equilibrium distribution, one obtains 
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FIG. 9: The perturbation \5f^ (bottom) is shown in a rotated frame such that f'^'^ is diagonal; the 
resulting oscillatory and damped behavior (top) of the element /12 of / in the box basis is due to 
the perturbation. The left and right figures correspond to the points marked (a) and (b) in Fig. ^ 
respectively. 

finally the linear response equations for the fundamental modes 

+ /:«[5/(^)]t + z/:«[5/(^)]t, (52) 

where we have defined a linear response operator through an appropriate centered difference 
limit 



^ •'^ x^o 2A ^ ^ 

We solve Eq. (^2]) as an eigenvalue problem. In general, the eigenvalues are complex, 
with frequency and damping rate given by the real and the imaginary parts respectively. 
The eigenvalues appear as complex conjugate pairs. The eigenmodes corresponding to non- 
zero eigenvalues are Hermitian conjugates of each other and are traceless with normalization 
given by 

Tr{5/W5/Wn = 1- (54) 

The physical linear response mode is given by the quadrature components Sf^ and Sf^. 
There also exists a zero mode that has nonvanishing trace. The damping rates are all 
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FIG. 10: The perturbation X5f^ (bottom) is shown in a rotated frame such that /^'^ is diagonal; the 
resulting damped behavior (top) of the element /12 of / in the box basis is due to the perturbation. 
The left and right figures correspond to the points marked (c) and (d) in Fig. respectively. 

negative, thus confirming the stabihty of the collective modes. In Fig. we plot the positive 
frequency eigenvalues. The dotted lines correspond to the difference frequencies of the 
Hamiltonian Hhf- 

It is interesting to see how these different modes evolve in real time. For this we use 
the equilibrium distribution obtained in the previous section and perturb it with one of the 
quadrature components 



/ 



(55) 



where A = 0.2 determines the smallness of the perturbation. In particular we will consider 
the modes labeled by (a), (b), (c) and (d) in Fig. |^. 

The real-time response is shown by plotting the off-diagonal matrix element /12 of the 
single particle density matrix in the box basis as shown in Figs. |^ and 0. In Fig. |TI|, we plot 
the change in the total energy AE = E{f) — E{f^^) as a function of time. In cases (a) and 
(b) we see that the AE oscillates about zero and eventually goes to zero. This is expected 
because such a perturbation tends to create coherences, resulting in an energy change by 
the amount of i^coh (coherence energy) that would eventually decay down to zero. Similar 
damped behavior is observed for case (d). Such an oscillatory damped behavior of the total 
energy could be attributed to the Markov approximation in the collision integral Eq. (pT|). 
On the other hand, perturbations of the kind (c) increase the total particle number by the 
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FIG. 11: The change in the total energy IS.E = E{f) — E{f^'^) as the system relaxes to its 
new equilibrium. In the top figure, the solid and dashed lines correspond to cases (a) and (b), 
respectively. Similarly in the bottom figure, the solid and dashed lines correspond to cases (c) and 
(d), respectively. 

amount 6N = Tr{A5/i,*^} and hence result in a finite change in total energy. 

VI. CONCLUSION 

A non-Markovian version of the quantum kinetic theory is derived using the prescription 
of a nonequilibrium statistical operator method as outlined in Ref. [^. This theory is 
shown to conserve energy in the ?7 — limit. Inclusion of quasi-particle damping results in a 
description beyond the Born approximation with energy conservation to 0^ order even in the 
Markov limit. To obtain collision integrals that involve quasi-particle damping and conserves 
energy precisely, one will have to calculate the T-matrix in the full collision operator keeping 
terms of all orders in the interaction. 

We applied the generalized second order kinetic theory to the nonhomogeneous dilute 
Bose gas confined in a spherical box to numerically study the full non-equilibrium evolution 
of the system towards equilibrium. The self-consistent distribution f^'^ thus obtained is very 
close to the Bose-Einstein distribution as shown in Fig. |[ We also observe a significant 
Hartree-Fock self energy shift which depends on the single particle distribution function /. 
The form of the damping function is important in determining the line shape. Particularly, 
the function with a 1/ cosh type of behavior is found to be appropriate and gives improved 
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energy conservation due to smaller initial correlation effect. 

The importance of such a real-time calculation is apparent from the full real-time response 
calculation, where we have calculated the damping rates and frequencies. These damping 
rates correspond to a shorter time scale compared to the equilibration time scale, which 
depend on rates in and out of the various levels. 

This simple model of a spherical trap can be easily extended to a more realistic situation 



of a harmonic trap. As in Refs. [T^, ^ |33[, the condensed component can be easily included 
by introducing a symmetry broken mean-field, = (a,), as one of the relevant observables 
and Hartree-Fock-Bogoliubov quasi-particle excitations. Such a calculation will allow us to 
make experimentally verifiable predictions of damping rates of collective excitations. One 
can also explore the possibility of including a time dependent potential or an external force 
term to selectively excite one or more of the collective modes. These possibilities will be 



explored in a forthcoming article [35 



Acknowledgment s 

We are extremely thankful to J. Cooper for the numerous and invaluable discussions. S.B. 
and M.H. acknowledge the financial support from the U. S. Department of Energy, Office 
of Basic Energy Sciences via the Chemical Sciences, Geosciences and Biosciences Division. 
R.W. gratefully receives financial support from the Austrian Academy of Sciences through 
an APART grant. 



APPENDIX A: REFERENCE DISTRIBUTION 

The reference distribution cr|° | of Eq. {^) is parameterized through its expectation values 
in Eq. ([Tl|). From the structure of this quantum Gaussian operator, it follows that 

<Mo)/iCt(,.o)} = UHtM^t^ymto)- (Al) 

In the above, U represents the single particle propagator of Eq. ( [TB| ) acting in many-particle 
Fock-space and /C is the corresponding single particle Hilbert-space propagator of Eq. (p3|). 
This condition implies that 

a,,a(°)Tr{[^W,7.]} = -[i^(°),a(°)] (A2) 
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and was used to obtain Eq. (p!7|). 



APPENDIX B: MATRIX ELEMENTS 



The matrix element of Eq. ( ^Of ) can be evaluated easily by expanding the sin functions into 
CO- and counter-propagating complex exponents and by an additional partial integration. 
This results in eight separate terms, i.e., 

■ AH dx 

(j)^^^^/as = — sin (i x) sin (j x) sin (k x) sin (I x) — , 
TT Jq x^ 

= F{t + j + k-l) + F{i + j-k + l) 
+ F{i - j + k + I) + F{i - j - k - I) 

- F{i + j - k - I) - F{i - J + k - I) 

- F{i-j-k + l)-F{i + j + k + l), (Bl) 

where 

F{n) = ^[cos (mi) +mi Siirnr)], (B2) 
Si{z) = f^^^dt. (B3) 

An asymptotic expansion of the sine integral leads to the following approximation that is 
correct at the 1% level, i.e., 

F(0) = ^, (B4) 

^, , \ ,TX , , sin (riTr) 2cos('n7r), 
F(n > 0) ^ — [-|n7r| ^ — + 
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